In [129]:
# Import libraries
import numpy as np
import pandas as pd; import geopandas as gpd
import libpysal as ps
# visualization, basemap, interactive mapping
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
from shapely.geometry import Point
import contextily as cx
import fiona
%matplotlib inline
import pydeck as pdk
import folium
import osmapi as osm
import osmnx as ox
from colorcet import fire
# For API
import requests
import json
import shapely
import pprint as pp
# pip install sodapy for downloading CalEnvScreen 3.0 data
from sodapy import Socrata
#For regression and GWR:
import statsmodels.api as sm
import numpy as np
from patsy import dmatrices
import esda as esda
import splot as splot
from splot.libpysal import plot_spatial_weights
from splot.esda import moran_scatterplot, plot_moran
from mgwr.gwr import GWR, MGWR
from mgwr.sel_bw import Sel_BW
from mgwr.utils import compare_surfaces, truncate_colormap
# ------- Data discription-------
# State-wide (CA): geographic, health status, environment and traffic, socio-economic status, ethnicity
# City-wide (SF): Physical activity (bike sharing, bike station, park, tree, sidewalk)
# Sources:
# https://oehha.ca.gov/calenviroscreen/report/calenviroscreen-30
# https://www2.census.gov/geo/tiger/GENZ2018/shp/
# https://data.sfgov.org/resource/7jbp-yzp3.json (bike sharing)
# https://opendata.arcgis.com/datasets/b8ad067022374acd966db3026ddab6fc_0.geojson (bike rack)
# https://data.sfgov.org/Transportation/Map-of-Sidewalk-Widths/ygcm-bt3x (sidewalk)
# ------- (Geo)Dataframe naming-------
# State level:
# df1: CA geographic data
# df2: CalEnviroScreen 3,0 data, including health outcome (e.g. asthma), environs (e.g. traffic, PM2.5, O3),
# socio-eco status (e.g. popverty, ethnicity), etc
# df3: merge df1 and df2
# City level:
# df5: sidewalk
# df6: bike
# df7: bike sharing
# ------ presentation format -------
# I referenced the following webpage to produce the slideshow for the recorded presentation
# https://medium.com/learning-machine-learning/present-your-data-science-projects-with-jupyter-slides-75f20735eb0f
Data import¶
In [2]:
# import geo data for CA and name it df1
df1 = gpd.read_file('data/cb_2018_06_tract_500k/cb_2018_06_tract_500k.shp')
In [3]:
# display first 3 columns of df1
df1.head(3)
Out[3]:
In [4]:
# plot the geo data
df1.plot(color='none', edgecolor='gray', linewidth=.2, figsize=(5,10))
Out[4]:
In [5]:
# check the coordinate system
df1.crs
Out[5]:
In [6]:
# import CalEnvScreen 3.0 Data and name it df2, by following the instructions provided by the data source
# https://www.opendatanetwork.com/dataset/www.transparentrichmond.org/wi9v-3g3n
client = Socrata("www.transparentrichmond.org", None)
# Returned as JSON from API
results = client.get("wi9v-3g3n", limit=20000)
# Convert to pandas DataFrame
df2 = pd.DataFrame.from_records(results)
In [7]:
# display data, including health outcome (e.g. asthma), environs (e.g. traffic, PM2.5, O3, diesel PM),
# socioethnic status (e.g. popverty, education level, ethnicity)
df2.head(2)
Out[7]:
In [8]:
# display all the columns (i.e. variables)
df2.columns
Out[8]:
The variables of interest to this project include geographic data ('tract' and 'zip') as well as, those variables within four main dimensions, health status ('asthma','lbw','cvd'), environment and traffic ('traffic','diesel,'pm'), socio-economic status ('edu','ling','pov') and ethnicity ('white_pct','african_am', 'asian_amer'). For each variable employed in this project, the description is listed as follows:
'zip':Postal ZIP Code that the census tract falls within
'asthma': Age-adjusted rate of emergency department visits for asthma.
'lbw': Percent low birth weight
'cvd': Percent low birth weight
'traffic': Traffic density, in vehicle-kilometers per hour per road length, within 150 meters of the census tract boundary
'diesel': Diesel PM emissions from on-road and non-road sources
'pm': Annual mean PM 2.5 concentrations
'edu': Percent of population over 25 with less than a high school education
'ling': Percent limited English speaking households
'pov': Percent of population living below two times the federal poverty level
'white_pct': Percent of white households
'african_am': Percent of african households
'asian_amer': Percent of asian households
The description for other variables can be found via the link: https://oehha.ca.gov/calenviroscreen/report/calenviroscreen-30.
In [9]:
# Extract variables of interest
df2=df2[['tract','zip','pm','traffic', 'diesel','asthma','lbw','cvd','edu','ling','pov','white_pct','african_am','asian_amer']]
In [10]:
# (see the following two lines)
# Found two already imported dataset use different naming of polygon identification
In [11]:
df1.GEOID[:5]
Out[11]:
In [12]:
df2.tract[:5]
Out[12]:
In [13]:
# For df2, remove the last two character (e.g. ".0") in the string
for i in range(len(df2['tract'])):
if df2['tract'].iloc[i].find('.') != -1:
df2['tract'].iloc[i] = df2['tract'].iloc[i][:-2]
df2['tract'][:5]
Out[13]:
In [14]:
# For df1, remove the first character (i.e. "0")
df1.GEOID=df1.GEOID.map(lambda x: x[1:])
df1['GEOID'][:5]
Out[14]:
In [15]:
# merge the aforementioned two dataset (i.e. df1 and df2) as a new dataframe called df3
df3=pd.merge(df1, df2, left_on='GEOID', right_on='tract')
df3.head(2)
Out[15]:
In [16]:
# (See the following two lines)
# make some thematic maps with columns 'traffic density' and 'households living with poverty'
In [17]:
df3["traffic"] = pd.to_numeric(df3["traffic"], downcast="float")
base = df3.plot(color='none', edgecolor='gray', linewidth=.2, figsize=(4,3))
df3.plot(ax=base, column='traffic',legend=True)
# The result is a bit odd, which would be solved afterwards
Out[17]:
In [18]:
df3["pov"] = pd.to_numeric(df3["pov"], downcast="float")
base = df3.plot(color='none', edgecolor='gray', linewidth=.2, figsize=(4,3))
df3.plot(ax=base, column='pov',legend=True)
# The result makes sense and I make a map for each variable in the next line of code and provide some descriptions.
Out[18]:
CA-level Chloropleth Maps¶
In [19]:
# Plot all variables with chloropleth maps for 12 variables within the four dimensions
# (i.e. health, socioeconomic, ethnicity, environs)
# cannot see a clear pattern of traffic density
f,(ax1,ax2,ax3) = plt.subplots(3,4,figsize=(30,16))
f.suptitle("Chloropleth Maps", fontsize=32)
# Health
df3["cvd"] = pd.to_numeric(df3["cvd"], downcast="float")
df3["lbw"] = pd.to_numeric(df3["lbw"], downcast="float")
df3["asthma"] = pd.to_numeric(df3["asthma"], downcast="float")
# Socio-economic
df3["pov"] = pd.to_numeric(df3["pov"], downcast="float")
df3["edu"] = pd.to_numeric(df3["edu"], downcast="float")
df3["ling"] = pd.to_numeric(df3["ling"], downcast="float")
# Ethnicity
df3["white_pct"] = pd.to_numeric(df3["white_pct"], downcast="float")
df3["asian_amer"] = pd.to_numeric(df3["asian_amer"], downcast="float")
df3["african_am"] = pd.to_numeric(df3["african_am"], downcast="float")
# Environs
df3["traffic"] = pd.to_numeric(df3["traffic"], downcast="float")
df3["pm"] = pd.to_numeric(df3["pm"], downcast="float")
df3["diesel"] = pd.to_numeric(df3["diesel"], downcast="float")
df3.plot(ax=ax1[0], column='cvd',legend=True, cmap="cet_CET_L8")
ax1[0].set_title("Cardiovascular Disease", fontsize=16)
df3.plot(ax=ax2[0], column='asthma',legend=True, cmap="cet_CET_L8")
ax2[0].set_title("Asthma", fontsize=16)
df3.plot(ax=ax3[0], column='lbw',legend=True, cmap="cet_CET_L8")
ax3[0].set_title("Low birth weight", fontsize=16)
df3.plot(ax=ax1[1], column='white_pct',legend=True, cmap='cet_fire')
ax1[1].set_title("% White Pop.", fontsize=16)
df3.plot(ax=ax2[1], column='asian_amer',legend=True, cmap='cet_fire')
ax2[1].set_title("% Asian American Pop.", fontsize=16)
df3.plot(ax=ax3[1], column='african_am',legend=True, cmap='cet_fire')
ax3[1].set_title("% African American Pop.", fontsize=16)
df3.plot(ax=ax1[2], column='pov',legend=True, cmap='cet_CET_D6')
ax1[2].set_title("poverty", fontsize=16)
df3.plot(ax=ax2[2], column='edu',legend=True, cmap='cet_CET_D6')
ax2[2].set_title("# below High School Edu", fontsize=16)
df3.plot(ax=ax3[2], column='ling',legend=True, cmap='cet_CET_D6')
ax3[2].set_title("Language barrier", fontsize=16)
df3.plot(ax=ax1[3], column='traffic',legend=True,cmap='cet_CET_D4')
ax1[3].set_title("traffic volume", fontsize=16)
df3.plot(ax=ax2[3], column='pm',legend=True, cmap='cet_CET_D4')
ax2[3].set_title("PM2.5", fontsize=16)
df3.plot(ax=ax3[3], column='diesel',legend=True, cmap='cet_CET_D4')
ax3[3].set_title("Diesel PM", fontsize=16);
Descriptive statistics (box plots and histograms)¶
In [20]:
f=plt.figure(figsize=(18, 6))
f.suptitle("Box Plots", fontsize=24)
plt.subplot(261)
plt.boxplot(df3.cvd)
plt.title('Cardiovascular disease', fontsize=10)
plt.subplot(262)
plt.boxplot(df3.lbw)
plt.title('Low birth weight', fontsize=10)
plt.subplot(263)
plt.boxplot(df3.asthma)
plt.title('Asthma', fontsize=10)
plt.subplot(264)
plt.boxplot(df3.pov)
plt.title('poverty', fontsize=10)
plt.subplot(265)
plt.boxplot(df3.edu)
plt.title('Education level', fontsize=10)
plt.subplot(266)
plt.boxplot(df3.ling)
plt.title('Linguistic barrier', fontsize=10)
plt.subplot(267)
plt.boxplot(df3.white_pct)
plt.title('% White households', fontsize=10)
plt.subplot(268)
plt.boxplot(df3.asian_amer)
plt.title('% Asian Am. households', fontsize=10)
plt.subplot(269)
plt.boxplot(df3.african_am)
plt.title('% African Am. households', fontsize=10)
plt.subplot(2,6,10)
plt.boxplot(df3.traffic)
plt.title('Traffic density', fontsize=10)
plt.subplot(2,6,11)
plt.boxplot(df3.diesel)
plt.title('Diesel PM', fontsize=10)
plt.subplot(2,6,12)
plt.boxplot(df3.pm)
plt.title('PM2.5 concentration', fontsize=10);
In [21]:
df3['traffic'].describe()
Out[21]:
In [22]:
# remove the outliers in 'traffic density'
low=df3['traffic'].quantile(.01)
high=df3['traffic'].quantile(.99)
high
Out[22]:
In [23]:
df3_f = df3[(df3['traffic'] < high) & (df3['traffic'] > low)]
df3_f['traffic'].describe()
Out[23]:
In [24]:
f=plt.figure(figsize=(18, 6))
f.suptitle("Box Plots", fontsize=24)
plt.subplot(261)
plt.boxplot(df3.cvd)
plt.title('Cardiovascular disease', fontsize=10)
plt.subplot(262)
plt.boxplot(df3.lbw)
plt.title('Low birth weight', fontsize=10)
plt.subplot(263)
plt.boxplot(df3.asthma)
plt.title('Asthma', fontsize=10)
plt.subplot(264)
plt.boxplot(df3.pov)
plt.title('poverty', fontsize=10)
plt.subplot(265)
plt.boxplot(df3.edu)
plt.title('Education level', fontsize=10)
plt.subplot(266)
plt.boxplot(df3.ling)
plt.title('Linguistic barrier', fontsize=10)
plt.subplot(267)
plt.boxplot(df3.white_pct)
plt.title('% White households', fontsize=10)
plt.subplot(268)
plt.boxplot(df3.asian_amer)
plt.title('% Asian Am. households', fontsize=10)
plt.subplot(269)
plt.boxplot(df3.african_am)
plt.title('% African Am. households', fontsize=10)
plt.subplot(2,6,10)
plt.boxplot(df3_f.traffic)
plt.title('Traffic density', fontsize=10)
plt.subplot(2,6,11)
plt.boxplot(df3.diesel)
plt.title('Diesel PM', fontsize=10)
plt.subplot(2,6,12)
plt.boxplot(df3.pm)
plt.title('PM2.5 concentration', fontsize=10);
plt.show();
In [25]:
# another way for descriptive statistics, i.e. histograms for observing data distribution
fig, ax = plt.subplots(2, 6,figsize=(20,5))
fig.suptitle("Histograms", fontsize=20)
m=12
# starting from 'pm' that is indexed as 12
for i in range(2):
for j in range(6):
df3_f.hist(column = df3_f.columns[m], ax=ax[i,j], figsize=(20, 18))
m+=1
;
Out[25]:
In [26]:
# We have three geo-levels, i.e. tract, zip, city (SF).
# For data reduction, I used the groupby operation to produce mean values of some variables for every zip code
df3_1=df3_f.groupby(df3_f['zip']).mean()
df3_1.head()
Out[26]:
In [27]:
len(df3)
Out[27]:
In [28]:
# Reduce data from 8034 to 1334 rows
len(df3_1)
Out[28]:
In [29]:
# white people has low lingustic barrier, as expected
lm=sns.lmplot("white_pct", "ling", data=df3_1, fit_reg=True, scatter=True, scatter_kws={"s": .1})
lm.fig.set_figwidth(3)
lm.fig.set_figheight(3)
In [30]:
# lower educational level positively correlate with low lingustic barrier, as expected
lm=sns.lmplot("edu", "ling", data=df3_1,scatter_kws={"s": .1})
lm.fig.set_figwidth(3)
lm.fig.set_figheight(3)
In [31]:
# diesel PM positively correlate with traffic that is one of the main sources of the diesel PM
lm=sns.lmplot("diesel", "traffic", data=df3_1,scatter_kws={"s": .1})
lm.fig.set_figwidth(3)
lm.fig.set_figheight(3)
In [32]:
# less affluent households tend to have health issue concerning low birth weight
lm=sns.lmplot("pov", "lbw", data=df3_1,scatter_kws={"s": 1})
lm.fig.set_figwidth(3)
lm.fig.set_figheight(3)
In [33]:
# places with more white households tend not to have asthmatic issue
lm=sns.lmplot("white_pct", "asthma", data=df3_1,scatter_kws={"s": 1})
lm.fig.set_figwidth(3)
lm.fig.set_figheight(3)
CA-level correlation analysis results¶
In [34]:
dfData = df3_1.iloc[:,2:].corr()
plt.subplots(figsize=(8, 8))
sns.heatmap(dfData, annot=True,vmax=1, vmin=-1, square=False, cmap="RdBu", alpha=.9, linewidths=2, linecolor='white', cbar=True)
plt.suptitle('Correlation Analysis', fontsize=20);
In [35]:
dfData = df3_1.iloc[:,5:11].corr()
plt.subplots(figsize=(6, 6))
mask = np.zeros_like(dfData)
mask[np.triu_indices_from(mask)] = True
sns.heatmap(dfData, annot=True,vmax=1, vmin=-1,mask=mask, square=False, cmap="RdBu", alpha=.9, linewidths=2, linecolor='white', cbar=True)
plt.suptitle('Correlation Analysis', fontsize=15);
In [36]:
# tracts with more white households tend not to bear social and health burdens
sns.pairplot(df3_1,vars=["edu","ling","pov",'asthma',"white_pct"],height=1,diag_kind="kde",corner=True,plot_kws=dict(marker="+", linewidth=.1))
Out[36]:
Data extraction¶
In [37]:
# Extract data with the label '075' in the column 'COUNTYFP' and display the first two rows.
SF = df3[df3['COUNTYFP']=='075']
SF.head(2)
Out[37]:
In [38]:
# There are 195 data points left.
len(SF)
Out[38]:
In [39]:
# Plot the geodataframe SF
SF.plot(color='none', edgecolor='gray', linewidth=.2, figsize=(5,10))
Out[39]:
In [40]:
# Make box plots again to make sure no erroneous result
# No outliers found
In [41]:
plt.figure(figsize=(16, 4))
plt.subplot(261)
plt.boxplot(df3.cvd)
plt.title('Cardiovascular disease', fontsize=10)
plt.subplot(262)
plt.boxplot(df3.lbw)
plt.title('Low birth weight', fontsize=10)
plt.subplot(263)
plt.boxplot(df3.asthma)
plt.title('Asthma', fontsize=10)
plt.subplot(264)
plt.boxplot(df3.pov)
plt.title('poverty', fontsize=10)
plt.subplot(265)
plt.boxplot(df3.edu)
plt.title('Education level', fontsize=10)
plt.subplot(266)
plt.boxplot(df3.ling)
plt.title('Linguistic barrier', fontsize=10)
plt.subplot(267)
plt.boxplot(df3.white_pct)
plt.title('% White households', fontsize=10)
plt.subplot(268)
plt.boxplot(df3.asian_amer)
plt.title('% Asian Am. households', fontsize=10)
plt.subplot(269)
plt.boxplot(df3.african_am)
plt.title('% African Am. households', fontsize=10)
plt.subplot(2,6,10)
plt.boxplot(df3_f.traffic)
plt.title('Traffic density', fontsize=10)
plt.subplot(2,6,11)
plt.boxplot(df3.diesel)
plt.title('Diesel vehicle', fontsize=10)
plt.subplot(2,6,12)
plt.boxplot(df3.pm)
plt.title('PM2.5 concentration', fontsize=10);
SF-level correlation analysis results¶
In [42]:
dfData = SF[['pm','traffic','diesel','asthma','lbw','cvd','edu','ling','pov','white_pct','african_am','asian_amer']].corr()
plt.subplots(figsize=(8,8))
sns.heatmap(dfData, annot=True, vmax=1, square=False, cmap="RdBu", alpha=.9, linewidths=2, linecolor='white', cbar=True)
plt.show()
In [43]:
dfData = SF[['edu','ling','pov','white_pct','asian_amer','african_am']].corr()
plt.subplots(figsize=(4, 4))
sns.heatmap(dfData, annot=True, vmax=1, square=False, cmap="RdBu", alpha=.9, linewidths=2, linecolor='white', cbar=True)
plt.show()
SF-level Chloropleth maps¶
In [44]:
f,(ax1,ax2,ax3) = plt.subplots(3,4,figsize=(30,16))
f.suptitle("Chloropleth Maps", fontsize=32)
# # Health
# SF["cvd"] = pd.to_numeric(SF["cvd"], downcast="float")
# SF["lbw"] = pd.to_numeric(SF["lbw"], downcast="float")
# SF["asthma"] = pd.to_numeric(SF["asthma"], downcast="float")
# # Socio-economic
# SF["pov"] = pd.to_numeric(SF["pov"], downcast="float")
# df3["edu"] = pd.to_numeric(SF["edu"], downcast="float")
# df3["ling"] = pd.to_numeric(SF["ling"], downcast="float")
# # Ethnicity
# df3["white_pct"] = pd.to_numeric(df3["white_pct"], downcast="float")
# df3["asian_amer"] = pd.to_numeric(df3["asian_amer"], downcast="float")
# df3["african_am"] = pd.to_numeric(df3["african_am"], downcast="float")
# # Environs
# df3["traffic"] = pd.to_numeric(df3["traffic"], downcast="float")
# df3["pm"] = pd.to_numeric(df3["pm"], downcast="float")
# df3["diesel"] = pd.to_numeric(df3["diesel"], downcast="float")
SF.plot(ax=ax1[0], column='cvd',legend=True, cmap="cet_fire")
ax1[0].set_title("Cardiovascular Disease", fontsize=16)
SF.plot(ax=ax2[0], column='asthma',legend=True, cmap="cet_fire")
ax2[0].set_title("Asthma", fontsize=16)
SF.plot(ax=ax3[0], column='lbw',legend=True, cmap="cet_fire")
ax3[0].set_title("Low birth weight", fontsize=16)
SF.plot(ax=ax1[1], column='pov',legend=True, cmap='cet_CET_L17_r')
ax1[1].set_title("poverty", fontsize=16)
SF.plot(ax=ax2[1], column='edu',legend=True, cmap='cet_CET_L17_r')
ax2[1].set_title("# below High School Edu", fontsize=16)
SF.plot(ax=ax3[1], column='ling',legend=True, cmap='cet_CET_L17_r')
ax3[1].set_title("Language barrier", fontsize=16)
SF.plot(ax=ax1[2], column='white_pct',legend=True, cmap='copper')
ax1[2].set_title("% White Pop.", fontsize=16)
SF.plot(ax=ax2[2], column='asian_amer',legend=True, cmap='copper')
ax2[2].set_title("% Asian American Pop.", fontsize=16)
SF.plot(ax=ax3[2], column='african_am',legend=True, cmap='copper')
ax3[2].set_title("% African American Pop.", fontsize=16)
SF.plot(ax=ax1[3], column='traffic',legend=True,cmap='cet_CET_D6_r')
ax1[3].set_title("traffic volume", fontsize=16)
SF.plot(ax=ax2[3], column='pm',legend=True, cmap='cet_CET_D6_r')
ax2[3].set_title("PM2.5", fontsize=16)
SF.plot(ax=ax3[3], column='diesel',legend=True, cmap='cet_CET_D6_r')
ax3[3].set_title("Diesel cars", fontsize=16);
Adding another dimension of variables concerning physical activity (walking/biking facilities)¶
1. Sidewalk¶
In [45]:
# load sidewalk data
df5 = pd.read_csv("data/MTA.sidewalk_widths.csv")
df5.head(2)
Out[45]:
In [46]:
# 16174 raws in total
len(df5)
Out[46]:
In [47]:
# the shape column is not really the "shape" for geopandas
type(df5['shape'][0])
Out[47]:
In [48]:
# try to transform the string type to shapely.geometry.linestring type
# print only one result for checking
shapely.wkt.loads(df5['shape'][0])
Out[48]:
In [49]:
# transforming the data type for all rows
df5['shape_'] = df5['shape'].apply(lambda x: shapely.wkt.loads(x))
In [50]:
df5.head(3)
Out[50]:
In [51]:
# transforming dataframe (df5) to GEOdataframe (gdf5)
gdf5 = gpd.GeoDataFrame(df5, crs=df1.crs, geometry=df5['shape_'])
gdf5.head(2)
Out[51]:
In [52]:
# displaying the sidewalk width range
gdf5['SIDEWALK_F'].describe()
Out[52]:
In [53]:
# filter out those less than 15m (below 75% perceptile)
base = SF.plot(color='none', edgecolor='r', linewidth=.1, figsize=(6,6))
gdf5[gdf5['SIDEWALK_F']>=15].plot(ax=base,edgecolor='gray', linewidth=.2, figsize=(5,5), legend=True)
Out[53]:
In [54]:
# intersect SF (main dataset) with gdf5 (sidewalk data) and create a new geodataframe called sidewalk
sidewalk=gpd.overlay(gdf5[gdf5['SIDEWALK_F']>=15], SF, how='intersection')
In [55]:
# now every sidwalk data point have a reference column for "groupby"ing back
sidewalk.head(1)
Out[55]:
In [56]:
sidewalk['sidewalk_length']=sidewalk.length
# calculating sidewalk length and store the values to a new column called sidewalk_length
In [57]:
# "groupby"ing back using "TRACTCE"
s=sidewalk.groupby('TRACTCE').sum()
In [58]:
len(SF)
Out[58]:
In [59]:
# there are 16 census tracts without any sidewalk, because there should be 195 rows
len(s)
Out[59]:
In [60]:
# Merge only the column 'sidewalk_length' back to SF
SF = SF.merge(s['sidewalk_length'],how='left', left_on='TRACTCE', right_index=True)
In [61]:
SF.plot(column='sidewalk_length',legend=True)
Out[61]:
In [62]:
# Show the null values in the 'sidewalk_length' column
SF[SF['sidewalk_length'].isnull()].head(2)
Out[62]:
In [63]:
# fill up the null values with zeros
SF['sidewalk_length'].fillna(0, inplace=True)
In [64]:
# check the null values have been filled with zeros
SF[SF['sidewalk_length'].isnull()]
Out[64]:
In [65]:
# display the data again
# found that sidewalk length is greater in western part of SF
SF.plot(column='sidewalk_length',legend=True)
Out[65]:
2. Bike rack¶
In [66]:
# load bike rack data
# looking at the json format for how I should extract data I need
url = "https://data.sfgov.org/resource/hn4j-6fx5.geojson"
response = requests.get(url)
results = response.text
data = json.loads(results)
data['features'][0]
Out[66]:
In [67]:
# Extract lat/lon data of bike racks
df6 = pd.DataFrame(columns=['lon','lat'])
for i in range(len(data['features'])):
df6.loc[i] = data['features'][i]['geometry']['coordinates']
df6.head
Out[67]:
In [68]:
# display the dataframe
shape = [Point(xy) for xy in zip(df6['lon'],df6['lat'])]
gdf6 = gpd.GeoDataFrame(df6, crs=df1.crs, geometry=shape)
gdf6.head()
Out[68]:
In [69]:
# plot the bike rack data
base = SF.plot(color='none', edgecolor='gray', linewidth=.1, figsize=(5,5))
gdf6.plot(ax=base, color='grey', markersize=.1)
Out[69]:
In [70]:
# using spatial join and groupby operation to prepare data
# counting bike rack numbers within each census tract
c=gpd.sjoin(SF, gdf6).groupby('GEOID').count()
In [71]:
# Randomly pick 'STATEFP' and store the count value to a new column called 'bikerack_count'
c['bikerack_count']=c['STATEFP']
c.head(1)
Out[71]:
In [72]:
# merge the 'bikerack_count' and SF main dataset
SF=SF.merge(c['bikerack_count'],how='left',on='GEOID')
In [73]:
# fill up null values with zero to represent no bike rack found
SF['bikerack_count'].fillna(0, inplace=True)
In [74]:
# found that bike racks concentrate on the northeasten part of SF
SF.plot(column='bikerack_count')
Out[74]:
3. Bike sharing station¶
In [75]:
# load bikesharing station data
url="https://data.sfgov.org/resource/7jbp-yzp3.json"
response = requests.get(url)
results = response.text
# data = json.loads(results)
data= pd.read_json(results)
In [76]:
# transform the data's latlon to float type
# create a geodataframe
fil = lambda x: float(x)
shape = [Point(xy) for xy in zip(data.station_longitude.map(fil), data.station_latitude.map(fil))]
gdf7 = gpd.GeoDataFrame(data, crs=df1.crs, geometry=shape)
gdf7.head(2)
Out[76]:
In [77]:
# Drop missing data
gdf7=gdf7[gdf7.station_latitude!=0]
In [78]:
# display the result
# found bikesharing station concentrate northeastern part of SF, similar to bike rack facilities
base= SF.plot(color='none', edgecolor='gray', linewidth=.1, figsize=(4,4))
gdf7[gdf7["station_id_domo"].str.contains('SF')].plot(ax=base, markersize=.5)
Out[78]:
In [79]:
# create a new column to store a value that represents the minimum value to the closest bike sharing station
SF['min_dist_to_bs'] = SF.geometry.apply(lambda g: gdf7[gdf7["station_id_domo"].str.contains('SF')].distance(g).min())
In [80]:
SF.head(2)
Out[80]:
In [81]:
# display the result, accssibility to bike sharing station
SF.plot(column='min_dist_to_bs',legend=True)
Out[81]:
Preparing data for GWR¶
In [82]:
SF.crs
Out[82]:
In [83]:
SF_3857 = SF.to_crs("EPSG:3857")
SF_3857.crs
Out[83]:
In [84]:
# create two columns 'x' and 'y' to store x,y coordinates using centroid function
SF_3857["x"]=SF_3857.centroid.map(lambda p: p.x)
SF_3857["y"]=SF_3857.centroid.map(lambda p: p.y)
SF_3857.head(2)
Out[84]:
In [85]:
fig, ax = plt.subplots(figsize = (6, 6))
SF_3857.plot(ax=ax, **{'edgecolor': 'black', 'facecolor': 'white'})
SF_3857.centroid.plot(ax = ax, c = 'black', markersize=10)
plt.show()
In [86]:
#Create choropleth quantile maps of each variable (including three new physical activity variables)
f,(ax1,ax2,ax3) = plt.subplots(3,5,figsize=(32,16))
# Health
SF_3857.plot(column='asthma', legend=True, ax=ax1[0],
scheme="quantiles", k=5, cmap='GnBu',edgecolor='black')
ax1[0].set_title("Asthma", fontsize=16)
SF_3857.plot(column='cvd', legend=True, ax=ax2[0],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax2[0].set_title("Cardiovascular disease", fontsize=16)
SF_3857.plot(column='lbw', legend=True, ax=ax3[0],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax3[0].set_title("Low birth weight", fontsize=16)
# Env
SF_3857.plot(column='diesel', legend=True, ax=ax1[1],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax1[1].set_title("Diesel cars", fontsize=16)
SF_3857.plot(column='traffic', legend=True, ax=ax2[1],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax2[1].set_title("Traffic density", fontsize=16)
SF_3857.plot(column='pm', legend=True, ax=ax3[1],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax3[1].set_title("PM2.5", fontsize=16)
# Socioeconomic
SF_3857.plot(column='edu', legend=True, ax=ax1[2],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax1[2].set_title("# people with education level below highschool", fontsize=16)
SF_3857.plot(column='ling', legend=True, ax=ax2[2],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax2[2].set_title("Linguistic barrier", fontsize=16)
SF_3857.plot(column='pov', legend=True, ax=ax3[2],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax3[2].set_title("Poverty", fontsize=16)
# Ethnicity
SF_3857.plot(column='white_pct', legend=True, ax=ax1[3],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax1[3].set_title("% white households", fontsize=16)
SF_3857.plot(column='asian_amer', legend=True, ax=ax2[3],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax2[3].set_title("% asian am households", fontsize=16)
SF_3857.plot(column='african_am', legend=True, ax=ax3[3],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax3[3].set_title("% african am households", fontsize=16)
# Physical activity
SF_3857.plot(column='sidewalk_length', legend=True, ax=ax1[4],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax1[4].set_title("Sidewalk length", fontsize=16)
SF_3857.plot(column='bikerack_count', legend=True, ax=ax2[4],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax2[4].set_title("# bike racks", fontsize=16)
SF_3857.plot(column='min_dist_to_bs', legend=True, ax=ax3[4],
scheme='quantiles', cmap='GnBu', k=5, edgecolor='black')
ax3[4].set_title("Accssibility to bikesharing station", fontsize=16)
for i in range(5):
ax1[i].axis('off')
ax2[i].axis('off')
ax3[i].axis('off')
cx.add_basemap(ax1[i], crs=SF_3857.crs.to_string(), source=cx.providers.Stamen.TonerLite)
cx.add_basemap(ax2[i], crs=SF_3857.crs.to_string(), source=cx.providers.Stamen.TonerLite)
cx.add_basemap(ax3[i], crs=SF_3857.crs.to_string(), source=cx.providers.Stamen.TonerLite)
Results of SF-level correlation analysis including physical activity¶
In [87]:
dfData = SF[['traffic','diesel','pm','asthma','lbw','cvd','sidewalk_length', 'bikerack_count', 'min_dist_to_bs']].corr()
plt.subplots(figsize=(6, 6))
sns.heatmap(dfData, annot=True, vmax=1, square=False, cmap="RdBu", alpha=.9, linewidths=2, linecolor='white', cbar=True)
plt.show()
Global regression analaysis: steps and major findings¶
In [88]:
# First model - 12 variables with no interaction term
y, X = dmatrices('asthma ~ pm + traffic + diesel +edu +ling + pov + white_pct + african_am + asian_amer +sidewalk_length + bikerack_count + min_dist_to_bs',
data=SF_3857, return_type='dataframe')
mod = sm.OLS(y, X)
res = mod.fit()
residuals = res.resid
predicted = res.fittedvalues
observed = y
print(res.summary())
In [89]:
SF_3857['bike_diesel_rack']=SF_3857['diesel']*SF_3857['bikerack_count']
In [90]:
SF_3857['bike_diesel_bshare']=SF_3857['diesel']*SF_3857['min_dist_to_bs']
In [91]:
# Second model - 14 variables with two interaction terms
y, X = dmatrices('asthma ~ pm + traffic + diesel +edu +ling + pov + white_pct + african_am + asian_amer +sidewalk_length + bikerack_count + min_dist_to_bs + bike_diesel_rack + bike_diesel_bshare',
data=SF_3857, return_type='dataframe')
mod = sm.OLS(y, X)
res = mod.fit()
residuals = res.resid
predicted = res.fittedvalues
observed = y
print(res.summary())
In [92]:
# Third model (result of backward stepwise regression)
y, X = dmatrices('asthma ~diesel + white_pct + african_am + asian_amer + bikerack_count + bike_diesel_rack',
data=SF_3857, return_type='dataframe')
mod = sm.OLS(y, X)
res = mod.fit()
residuals = res.resid
predicted = res.fittedvalues
observed = y
print(res.summary())
Spatial autocorrelation and GWR analysis results:¶
Spatial autocorrelation¶
In [93]:
SF_3857.crs
SF_3857.reset_index(drop=True, inplace=True)
SF_3857.head(2)
# resetting value for caculating kernal weight
Out[93]:
In [94]:
w = ps.weights.Kernel.from_dataframe(SF_3857, k=5, fixed=False, function='gaussian')
plot_spatial_weights(w, SF_3857, figsize=(5,5));
# w = ps.weights.Queen.from_dataframe(SF_3857)
# plot_spatial_weights(w, SF_3857, figsize=(5,5));
In [95]:
# Global Moran's I
y, X = dmatrices('asthma ~ pm + traffic + diesel +edu +ling + pov + white_pct + african_am + asian_amer +sidewalk_length + bikerack_count + min_dist_to_bs',
data=SF_3857, return_type='dataframe')
moran = esda.Moran(y, w)
moran.I
Out[95]:
In [96]:
# Plotting the result
# A Moran's I of .82 is quite high
fig, ax = moran_scatterplot(moran, aspect_equal=True)
plt.show()
In [97]:
# calculating Moran_Local and plotting the result
moran_loc = esda.moran.Moran_Local(y, w)
fig, ax = splot.esda.moran_scatterplot(moran_loc, p=0.05)
ax.set_xlabel('Asthma')
ax.set_ylabel('Spatial Lag of Asthma')
plt.show()
In [98]:
# plotting the result on the map
# two hotspots and one large coldspot can be found, which are very similar to the residential patterns of white and Asian households
moran_loc = esda.moran.Moran_Local(y,w)
# splot.esda.lisa_cluster(moran_loc, SF_3857, p=0.05, alpha=.9, figsize = (5,5))
splot.esda.plot_local_autocorrelation(moran_loc, SF_3857, 'asthma', p=0.01, quadrant=1)
plt.show()
GWR¶
In [99]:
#Prepare dataset inputs for GWR
g_y = SF_3857['asthma'].values.reshape((-1, 1))
g_X = SF_3857[['diesel', 'white_pct', 'african_am', 'asian_amer', 'bikerack_count','bike_diesel_rack']].values
u = SF_3857['x']
v = SF_3857['y']
g_coords = list(zip(u, v))
In [100]:
#Calibrate a GWR model for Chicago dataset using computationally selected~bandwidth (optimizes using AIC)
gwr_selector = Sel_BW(g_coords, g_y, g_X)
gwr_bw = gwr_selector.search()
print('Computationally selected bandwidth: ',gwr_bw)
gwr_model = GWR(g_coords, g_y, g_X, gwr_bw)
gwr_results = gwr_model.fit()
print('Residual sum of squares: ', gwr_results.resid_ss)
In [101]:
gwr_results.summary();
# X1:'diesel', X2:'white_pct', X3:'african_am', X4:'asian_amer', X5:'bikerack_count', X6:'bike_diesel_rack'
In [102]:
#Prepare GWR results for mapping
#Add GWR parameters to GeoDataframe
SF_3857['gwr_intercept'] = gwr_results.params[:,0]
SF_3857['gwe_diesel'] = gwr_results.params[:,1]
SF_3857['gwr_white'] = gwr_results.params[:,2]
SF_3857['gwr_african'] = gwr_results.params[:,3]
SF_3857['gwr_asian'] = gwr_results.params[:,4]
SF_3857['gwr_bikerack_count'] = gwr_results.params[:,5]
#Obtain t-vals filtered based on multiple testing correction
gwr_filtered_t = gwr_results.filter_tvals()
In [103]:
# Plot estimators
fig, (axes1, axes2) = plt.subplots(nrows=2, ncols=3, figsize=(12,6), constrained_layout=True)
ax0 = axes1[0]
ax0.set_title('GWR Intercept Surface (BW: ' + str(gwr_bw) +')', fontsize=10)
ax1 = axes1[1]
ax1.set_title('GWR diesel Surface (BW: ' + str(gwr_bw) +')', fontsize=10)
ax2 = axes1[2]
ax2.set_title('GWR white Surface (BW: ' + str(gwr_bw) +')', fontsize=10)
ax3 = axes2[0]
ax3.set_title('GWR african Surface (BW: ' + str(gwr_bw) +')', fontsize=10)
ax4 = axes2[1]
ax4.set_title('GWR asian Surface (BW: ' + str(gwr_bw) +')', fontsize=10)
ax5 = axes2[2]
ax5.set_title('GWR bike rack count (BW: ' + str(gwr_bw) +')', fontsize=10);
SF_3857.plot('gwr_intercept', legend=True, cmap='seismic', ax=ax0, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,0] == 0).any():
SF_3857[gwr_filtered_t[:,0] == 0].plot(color='lightgrey', ax=ax0, **{'edgecolor':'black'})
SF_3857.plot('gwe_diesel', legend=True,cmap='seismic', ax=ax1, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,1] == 0).any():
SF_3857[gwr_filtered_t[:,1] == 0].plot(color='lightgrey', ax=ax1, **{'edgecolor':'black'})
SF_3857.plot('gwr_white', legend=True, cmap='seismic', ax=ax2, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,2] == 0).any():
SF_3857[gwr_filtered_t[:,2] == 0].plot(color='lightgrey', ax=ax2, **{'edgecolor':'black'})
SF_3857.plot('gwr_african', legend=True, cmap='seismic', ax=ax3, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,3] == 0).any():
SF_3857[gwr_filtered_t[:,3] == 0].plot(color='lightgrey', ax=ax3, **{'edgecolor':'black'})
SF_3857.plot('gwr_asian', legend=True, cmap='seismic', ax=ax4, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,4] == 0).any():
SF_3857[gwr_filtered_t[:,4] == 0].plot(color='lightgrey', ax=ax4, **{'edgecolor':'black'})
SF_3857.plot('gwr_bikerack_count', legend=True, cmap='seismic', ax=ax5, **{'edgecolor':'black', 'alpha':.65},legend_kwds={'shrink': .9})
if (gwr_filtered_t[:,5] == 0).any():
SF_3857[gwr_filtered_t[:,5] == 0].plot(color='lightgrey', ax=ax5, **{'edgecolor':'black'});
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
ax2.get_xaxis().set_visible(False)
ax2.get_yaxis().set_visible(False)
ax3.get_xaxis().set_visible(False)
ax3.get_yaxis().set_visible(False)
ax4.get_xaxis().set_visible(False)
ax4.get_yaxis().set_visible(False)
ax5.get_xaxis().set_visible(False)
ax5.get_yaxis().set_visible(False);
In [104]:
#Local model fit
SF_3857['R2'] = gwr_results.localR2
SF_3857.plot('R2', legend = True, figsize=(5,5),legend_kwds={'shrink': .8})
ax = plt.gca()
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
#plt.savefig('local_R2')
plt.show()
In [105]:
VIF = gwr_results.local_collinearity()[1]
VDP = gwr_results.local_collinearity()[3]
names = ['Diesel','White','African','Asian','Bike','interaction']
fig, ax = plt.subplots(1, 6, figsize = (20, 5))
for col in range(6):
SF_3857['vif'] = VIF[:, col]
SF_3857.plot('vif', ax = ax[col], legend = True, legend_kwds={'shrink': .5})
ax[col].set_title('VIF: ' + names[col])
ax[col].get_xaxis().set_visible(False)
ax[col].get_yaxis().set_visible(False)
fig, ax = plt.subplots(1, 6, figsize = (20, 5))
for col in range(6):
SF_3857['vdp'] = VDP[:, col]
SF_3857.plot('vdp', ax = ax[col], legend = True, legend_kwds={'shrink': .5})
ax[col].set_title('VDP: ' + names[col])
ax[col].get_xaxis().set_visible(False)
ax[col].get_yaxis().set_visible(False)
In [106]:
# Display bike rack data using Gridlayer from Mapbox
bike_rack = pdk.Layer(
'GridLayer',
df6,
get_position=['lon', 'lat'],
auto_highlight=True,
elevation_scale=150,
pickable=True,
elevation_range=[0, 30],
extruded=True,
opacity=.3,
coverage=.7)
view_state = pdk.ViewState(
longitude=-122.415,
latitude=37.8023,
zoom=9,
min_zoom=1,
max_zoom=15,
pitch=50,
bearing=30)
r = pdk.Deck(layers=[bike_rack], initial_view_state=view_state)
r.to_html('Grid.html')
Out[106]:
In [107]:
# Display bike sharing station data using Heatmaplayer from Mapbox
view = pdk.ViewState(
longitude=-122.455,
latitude=37.7423,
zoom=11,
min_zoom=1,
max_zoom=15)
COLOR_SCALE = [
[90, 249, 232],
[70, 204, 197],
[55, 150, 150],
[20, 100, 100],
]
bike_sharing = pdk.Layer(
"HeatmapLayer",
data=gdf7,
opacity=0.9,
get_position=["station_longitude", "station_latitude"],
aggregation='"MEAN"',
color_range=COLOR_SCALE,
threshold=1,
get_weight="dock_count",
pickable=True,
)
r = pdk.Deck(
layers=bike_sharing,
initial_view_state=view,
tooltip={"text": "Concentration of bike-sharing decks in orange"},
)
r.to_html("heatmap_layer.html")
Out[107]:
In [108]:
# Create a osmnx-graph called G_bike to represent bike network and displaying it
G_bike = ox.graph.graph_from_bbox(37.86,37.7,-122.325,-122.55,network_type='bike')
ox.plot.plot_graph(G_bike,ax=None, figsize=(8, 8), bgcolor='#111111',
node_color='w', node_size=.5, node_alpha=None, node_edgecolor='none', node_zorder=1,
edge_color='#999999', edge_linewidth=.1, edge_alpha=None, bbox=None)
Out[108]:
In [109]:
# input data for origin and destination
# here I assume that I am going to "Taiwan Restaurant, SF" from "Hunter's Point, SF"
origin_input = input()
destination_input = input()
In [110]:
# using geocoder to transform the string-type address to latlon type coordinates
origin = ox.geocoder.geocode(origin_input)
destination = ox.geocoder.geocode(destination_input)
In [111]:
# showing origin's latlon
origin
Out[111]:
In [112]:
# showing destination's latlon
destination
Out[112]:
In [113]:
# calculting the distance between the OD pair
ox.distance.euclidean_dist_vec(origin[0], origin[1], destination[0], destination[1])
Out[113]:
In [114]:
# getting the nearest edge ID for origin and destination
ii_bike=ox.distance.get_nearest_edge(G_bike, origin)
jj_bike=ox.distance.get_nearest_edge(G_bike, destination)
In [115]:
# showing edge ID
ii_bike
Out[115]:
In [116]:
# showing edge ID
jj_bike
Out[116]:
In [117]:
# showing edge IDs along the shortest path for the OD pair
path_bike = ox.distance.shortest_path(G_bike, ii_bike[0], jj_bike[0], weight='length')
path_bike[:10]
Out[117]:
In [118]:
# the total count of the edges on the shortest path
len(path_bike)
Out[118]:
In [119]:
# displaying all the info for all edges
edges = ox.utils_graph.get_route_edge_attributes(G_bike, path_bike, attribute=None, minimize_key='length', retrieve_default=None)
edges[:2]
Out[119]:
In [120]:
# Mapping the shortest route and diesel PM onto the folium platform
diesel_map = folium.Map([37.7556, -122.4399], zoom_start = 12)
tracts = SF.set_index('TRACTCE')
SF_wgs84 = SF.to_crs({'init': 'epsg:4326'})
folium.Choropleth(
geo_data=SF,
name='choropleth',
data=SF,
columns=['TRACTCE','diesel'],
key_on = 'feature.properties.{}'.format('TRACTCE'),
fill_color='YlOrBr',
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Diesel Emission'
).add_to(diesel_map)
route_diesel_map=ox.folium.plot_route_folium(G_bike, path_bike, route_map=diesel_map, fit_bounds=False,
popup_attribute=None, route_color='#5e8277', route_width=5, route_opacity=.9)
route_diesel_map
Out[120]:
In [121]:
# using nodeget function to derive the latlon for each edge and store them in a list called node_lst
# i.e. transforming edge ID "back" to latlon format
api = osm.OsmApi()
node_lst=[]
for i in range(len(path_bike)):
node_lst.append(api.NodeGet(path_bike[i]))
node_lst[:2]
Out[121]:
In [122]:
# create a list called node_latlon to store latlon data extracted from node_lst
node_latlon=[]
for i in range(len(node_lst)):
node_latlon.append(list([node_lst[i]["lon"], node_lst[i]["lat"]]))
In [123]:
type(node_latlon)
Out[123]:
In [124]:
# check the data format is correct as expected
node_latlon[0]
Out[124]:
In [125]:
# create a dataframe called d, with two columns which will be filled with coordiniate and timestamp data
# transform node_latlon to pd.Series
d = pd.DataFrame(columns=['coordinates','timestamps'])
pd.Series(node_latlon)
Out[125]:
In [126]:
# randomly insert 1 in the first row of d
# inserting latlon data into the column 'coordinates' of d
d.loc[0]=1
d['coordinates'][0]=node_latlon
d['coordinates']
Out[126]:
In [127]:
# assuming that each travel step takes 1 sec
# here, I did not use the "actual" estimate of travel time as I was not able to figure out how,
# so I just used 1, 2, 3 etc to represent the travel timestampstime_lst=[]
time_lst=[]
for i in range(len(node_latlon)):
time_lst.append(i)
d['timestamps'][0]=time_lst
d
Out[127]:
In [128]:
# mapping the shortest route onto another platform, Mapbox
layer = pdk.Layer(
"TripsLayer",
d,
get_path="coordinates",
get_timestamps="timestamps",
get_color=[200, 50, 150],
opacity=1,
width_min_pixels=10,
rounded=False,
trail_length=500,
current_time=300,
)
view_state = pdk.ViewState(latitude=37.7749295, longitude=-122.4194155, zoom=10, bearing=50, pitch=45)
r = pdk.Deck(layers=[layer,bike_rack,bike_sharing], initial_view_state=view_state)
# r = pdk.Deck(layers=layer, initial_view_state=view_state)
r.to_html("trips_layer.html")
Out[128]:
Thanks!